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Abstract 

Non-Markovian quantum state diffusion (NMQSD) is an exact method for calculating the reduced 
density matrix of an arbitrary subsystem interacting linearly with the radiation field. Applications 
of the theory have however been few due to the intractable nature of the variational-differential 
NMQSD evolution equation. Recently, we argued that the variational-differential equation can 
be rewritten as an integrodifferential equation which can be readily solved numerically. This 
manuscript provides an explicit derivation of the modified equations. Applications to intermittent 
fluorescence in 24 Mg + are discussed in detail. Earlier speculations that quantum jumps occur on all 
time scales are verified on a picosecond timescale. We show that a plot of the probability density of 
the signal vs signal strength shows the two characteristic peaks associated with the bright and dark 
manifolds, and that the ratio of the areas under the peaks is 16 as observed experimentally. We 
also show that the shape of this distribution is sensitive to bath memory, but has a mathematical 
form common to both the Markovian and non-Markovian cases. 
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I. INTRODUCTION 

Like quantum state diffusion^], non-Markovian quantum state diffusion (NMQSD) 0, 
was originally introduced as an exact computational method for finding the reduced 
density matrix of an arbitrary subsystem interacting linearly with the radiation field. In 
practice the evolution equation proved impossible to solve except in a few cases where exact 
solutions were already known 0]. Recently we reformulated NMQSD in terms of a solvable 
integrodifferential equation and demonstrated the use of the modified method by solving 
a number of example problems^. Here we present a detailed and general derivation of 
the modified equations for an arbitrary number of coupling operators. We also apply the 
resulting theory to intermittent fluorescence in driven 24 Mg + which arises due to quantum 
jumps between bright and dark electronic states. 



Ion trap experiments have become an important subfield of quantum optics. Theory 



■ • ■ M 

|8 1 and experiments |9l. Il0| on quantum jump phenomena in single ions are motivating further 
attempts to understand the measurement process HQ. 

Such ions have been used as 

models of quantum computers Q], and there are many interesting possibilities when optical 
lattices are emp loved [141] . In many cases the subsystem of interest in these experiments can 
be modeled as a driven few level system interacting linearly with the radiation field. NMQSD 
is an exact dynamical theory for such systems which should in principle make theoretical 
analysis of such experiments a simple task. Unfortunately, standard formulations of NMQSD 
rely on a stochastic variational-differential equation (VDE) for dynamical evolution which 
cannot be solved - even numerically - outside of a few special cases where solutions were 
already known. Some approximation schemes have recently been explored jisl. Il^ . but exact 
numerical schemes remain an important goal. 

In a recent manuscript we argued that NMQSD can be reformulated in terms of a solvable 
stochastic integrodifferential equation ji^]. We solved a number of example problems - some 
previously unsolvable - with the reformulated equations and found that accurate solutions 
to few level problems were readily obtainable. The efficiency of the method is partly a con- 
sequence of the availability of improved methods for solving stochastic differential equations 
(SDEs)0,Q- In this manuscript we apply the reformulated theory to 24 Mg + . Our results 



confirm speculation that quantum jumps occur on short time scales 



101 ] . We show that both 



Markovian and non-Markovian quantum state diffusion predict the two peaked distribution 
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seen in experiments. In addition both theories give the correct ratio for the areas under 
the peaks. We also find that the mathematical lineshape functions are very similar for both 
theories even though the peaks appear qualitatively different in the two cases. 

The derivation of the reformulated equations for NMQSD is outlined in section II. The 
equations and their numerical implementation are discussed in section III. In section IV we 
construct a three level model for 24 Mg + . Numerical results for application of NMQSD to 
the 24 Mg + model are discussed in section V. 

II. DERIVATION OF DYNAMICAL EQUATION 
A. Hermitian coupling 

For simplicity we will consider the simple special case where the total Hamiltonian is 

m m 

H tot = H + xJ29j(a] + aj) + J2 f^jOjCLj 
3=1 3=1 

= H + H C (1) 

and comment on the modifications necessary for more general Hamiltonians. We will also 
assume a bath temperature of K and hence an initial state of the form 

l*Uo)) = hM®|o)...®|o) (2) 

where |0) denotes the lowest eigenstate of a^hj. We will first show that the reduced density 

pit) = Tr 6at ,{e-^-*/ ?i |* tot (0))(* tot (0)|e^/ ?i } (3) 
can be rewritten as an average over diadics 

p{t) = M[\^ t )^ t \}. (4) 
To do this we will express the trace over bath modes as integrals over coherent states \(Xj), 

%K> = ( 5 ) 
for each oscillator. We do this by inserting closure relations 

/ rf2 «jK')(«il = h ( 6 ) 
3 



in Eq. (J3J), where cPctj = dRectj dlmaj/n. If x, y denote eigenvalues of x, and inserting 
closure relations J^° 00 dx'\x'){x'\ = 1 and J^° 00 dy'\y'){y'\ = 1, then it follows that matrix 
elements of the reduced density can be expressed as 



(x\p(t)\y) = dx' dy'i/) (x')il)o(y') / d 2 a x ... / d 2 a m 

J -co J— oo J J 

(x, ax,... , aJe-^/V, 0, . . . , 0) (y', 0, . . . , 0|e^ /a |y, a m >. (7) 

We now define dt = t/N and use the Trotter product formula 

e -iHtott/h _ jj m e -iHdt/h e -iH c dt/h e ~iHdt/h e ~iH c dt/h ^ 
N—+00 

where there are now 2N factors inside the limit. Note that for more general Hamiltonians 
where there are many coupling operators Lk which couple to different modes of the bath it 
is necessary to use a appropriate generalization of the Trotter product formula to separate 
the subsystem Hamiltonian H and coupling operators in separate factors. After this each 
coupling operator can be treated separately using techniques similar to those which follow. 

Inserting N — 1 closure relations for x in Eq. (JHJ) and inserting the result into the matrix 
element (x, ai, . . . , a m \e~ tHtott / h \x\ 0, . . . , 0) then gives 



/OO />00 
dxi . . . / dx 
-co J-oo 



in 



J~J ^ e -i[x N -ig](a.'l+d ] )+hujjd^d J ]dt/h e -i[x 1 g j (d' t j +d J )+huj J d^d : j]dt/h |g^ ^ 

The factors in the coherent state matrix element can now be combined using 

e Adt e Bdt _ e (A+B)dt+0{dt 2 ) 

and neglecting the 0(dt 2 ) terms. Introducing the usual path integral notation then gives 
(x, a m \e- iHt ° tt/h \x', 0, . . . , 0) = / V[x]e iS[x]/h J] (a, |e~' [ Jo &*t>9W^)+toj*frW\ ) 

Jx ' 3=1 

(11) 

where S[x) is the usual action. 

The coherent state matrix elements can now be found by considering their dynamics. 
Defining 



and using the facts that (aj|a] = a*(aj\ and (otj\dj = (<x,/2 + d/da*)(ctj\ it then follows 
that ipj(t,aj,a*) satisfies 

dipj(t,aj,a*)/dt = — (i/h){x t gj(a*+aj/2+d/da*)+hujja*(aj/2+d/da*)}'ipj{t,aj,a*) (13) 

with initial condition if)j(0, <x,-, a*) = (<x,|0) = e~ aja ^ 2 . Guessing a solution of the form 

ipj(t, otj, a*) = exp{— aja*/2 + Cj(i)a* + dj(t)} (14) 

and substituting into Eq. (|T3*j) one finds first order equations for unknowns Cj(t) and dj(t) 
which can be solved with initial conditions Cj(0) = and dj(0) = 0. The solutions are 

c At) = -H/%) f dt' Xtfgje-*"*-^ 
Jo 

dj(t) = -(1/h 2 ) f dt' f dt" xwg^e-^iV-*') (15) 



m 

aj a*/2 



and hence we may rewrite the matrix element in the form 

(x, a m \e- iA ^ n \x', 0, . . . , 0) = f T>[x]e iS ^ h f[ e 

Jx ' j=i 

e -(i/R) f* dt' It , Sj e-^( 4 -'\« e -(l/f» 2 ) f dt' f*' dt" x t ,x tll g]e- iu _ ^ 

A similar formula can be found for (y', 0, . . . , 0\e lHtott ^ h \y, cti, . . . , a m ) and when both are 
substituted into Eq. (J7J) we obtain 

/oo roc rx ry 

dx' / dy'Mx')%(y') / V[x] / 2%]e*H/» c -*%]/» 
-oo J— oo J x 1 J y f 

ft / d 2 aj e - aja *ie- {i/n) fo dt ' xt>9je-^ {t -^ a * ei i/ h) ft d t' vweW-^as 
■-1 J 

m 2 )f*dt df x t ,x t «^e- fa 'i(*'-*") e -(l/R a )/ t df eft" e^*'-'") _ ^ 



i=l 



The integrals over the real and imaginary parts of the a, are now just Gaussian integrals 
which can be performed analytically giving 



(x\p(t)\y)= dx' dy'Mx')r (y') V[x] V[y]e iS W n e- iS W h 

J —oo J —oo J x' Jy f 

e Io dt ' Jo dt " x t'Vt»**(t't") e - Jo Jo' dt" x t ,a V ,a(t',i») e - f o dt' J*' dt" y t ,y t „a*(t>,t») ^ 



where a(t,t') = (1/h 2 ) E™i {t ~ t ' ) - 

The generalization to non-zero temperatures requires consideration of all 
(x, aii, • • • , a m \e~ zHtott / h \x r , ni, . . . , n m ) matrix elements (and their y-y' counterparts) 



for all ni,...,n m quantum numbers. The equation (JT3)l for ipj(t,aj,a*) is unaltered but 
the initial condition is now A — e~ aja ^ 2 and so the correct ansatz is 

( a * _|_ b-(t)) nj 

ifjj(t,aj,a*) = — - — j= exp{— a.ja*/2 + Cj(t)a* + dj(t)} (19) 



and the solutions for bj(t), Cj(t) and dj(t) are 

bj(t) = -{i/h) f dt' xt.g^ 1 -^ 
c -(t) = -(i/h) f dt' xtgje-^'-V 



d 3 (t) = -(1/h 2 ) f dt' f dt" XfXtngle-^'-^ -irijUjt. (20) 

^ J 



The correct thermal weight for mode j is e ^"jAb^i _ exp(— htOj/kBT)) and introducing 
the notation Wj = exp^—huj/ksT) we get 



dx' / dy'Mx')r (y') / V[x] / P^e^NAe-^M/* 

-oo J —oo J x' J y 1 

i=i 

g -(lA 2 ) /„' df dt" x^x^e-^i <*'"*"> e -(l/^ 2 ) /„ «' Jo' *" W^jje^f-*") 

£ T-ifa + {i/h) dt'ytgje-^-t^-ii/h) dt'x.g^^ ))}"*. (21) 

Performing the sum explicitly then gives 

-oo J —oo Jx' Jy' 

3=1 J 

e -(l/h 2 ) ft dt' ft' dt" x^ng^ e -(l/ft 2 ) /„' dt' dt" ^fffe^*'-*") 

^(aj+tt/R) f^v^e-^^-^Xa^-ii/h) ft dt'x^g^i^) ^ 

or with some rearrangement 



J —oo J —oo J x' J y 1 



m „ 



-(i/h) ftdt' x^gje-^-^+ii/h) ft Of y^e-^-^Wjja* 



r (l/n 2 ) f* dt' dt" xwrfe-toS e ~m 2 ) jl dt> dt" wvtvgjeW-*') 

,(l/h 2 ) W] f* dt' f* dt"x t ,y t „g^ t '- t " ) _ (33) 



Again the integrals over real and imaginary parts of acj are Gaussian and can be done explic- 



itly. The result is again Eq. ()18j) where now a(t,t') = (l/h 2 )J2T=i9^[ co ^(^ J f) COSUJ j(t 



t') — i sin Uj(t — £')]. When the coupling operator is non-Hermitian both (JTSj) and the memory 
function must be modified^. 

Equation (|18J) was first obtained by Feynman and Vernon The first exponential factor 
in this expression couples x and y which means the path integrals are coupled and must be 
performed simultaneously. This coupling can be eliminated at the expense of introducing a 
complex stochastic process z t . Specifically, we use the identity 

e Io dt> Jo dt " x t'y t >>a*(t',t") _ M[eio dt ' ^t'^'+yf^)^ (24) 
first employed by Strunz[2j, where the mean over realizations of the noise is Gaussian 

M[F[z)) = J V[z] F[z) Ne~ 1°° dt ' Io°° dt " (25) 

and (3{t'X) is the functional inverse of a(t,t') (i.e. ^ dt' a(t,t')/3(t',t") = 5(t - t")). 
Equation ({21} can be proved using (|23|) by completing the square in the exponent. 
Finally, inserting ()24|) into Eq. ()18|) we can rearrange terms so that 

(x\m\y) = n(^t)(A\y)} (26) 

where 

( x \ij, t ) = [°° dx'M%') [ X V[x]e iS[x]/h efo dt ' x f z f e ~ lo dt ' Jo dt " (27) 

J — oo J x' 

To obtain a wave equation we differentiate with respect to time 

/OO /* 30 

dx'Mx') / V[x}e iS ^ h 
-oo J x' 

r f* dt' x f ,z,, - f* dt' f* dt" x f ,x.r,a(t' ,t") 



JO 
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and rewrite each term in terms of ipt- Noting that xt = x the second term presents no 
difficulties but the third does because of the delayed factor Xf. Non-Markovian quantum 
state diffusion resolves this problem by introducing a variational derivative via the identity 

5z t n 1 ; 

Thus, the linear NMQSD wave equation takes the form 

^ = -(i/h)Hip t + xip t z t - x f dt"a(t,t")^-. (30) 
at Jo ozt" 

Unfortunately, numerical methods for variational-differential equations like (J3*U|) have 

not yet been developed. This means that few problems can be solved using (J3U)) and its 

generalizations Q]. Recently, we have shown that there is a alternative way to formulate 

NMQSD in terms of a solvable integrodifferential equation^]. The key idea is that Eq. (|27|) 

can be rewritten in terms of a dynamical semigroup U(t, 0) via (x\ip t ) — ( x \U(t,0)\ipo) . To 

see that U (t, 0) does represent a semigroup note that Eq. (|27ji can be rewritten as 

J — OO J — OO JXi 

Xl V[x]e iS[x]/h efo dt ' x t' z t' e ~-fo l dt ' Jo dt " W'*^.*") (31) 



for any intermediate time tj. This means U(t,0) = U(t,ti)U(ti,Q) and hence we have a 
semigroup. Applying this to the second term in (|28j) enables us to write 



dU(t, 0)/dt = -(i/h)HU(t, 0) - ixU(t, 0)z t - x f dt"a(t, t")U(t, t")xU(t", 0). (32) 

Jo 

All considerations above remain unchanged under t — > — t and hence U(t,0)~ l exists. Using 
U{t,t") = E)'(t,0)£)'~ 1 (t / ',0) and changing notation to U t = U(t,0) we finally obtain 

dUJdt = -(i/h)HU t + xU t z t - xU t f dt"a(t, f)U^xU v > (33) 

Jo 

which is a closed integrodiffential equation for the propagator U t . 

As in the case of NMQSD we must now find a norm-preserving version of the theory. 
This is accomplished using the Girsonov transformation discussed at length in Ref . Q . The 
result is the evolution equation 

rt 
10 



dU t /dt = -(i/h)HU t - i(x - (x)t)U t (z t + / dt"a*(t,t")(x} t „) + C t U t 



-(x - (x)t)U t f dt"a(t, f)U^xU v > (34) 
Jo 



where the function Cj is given by 



C t = (^ \U}(x - (x) t )U t / dt"a(t,t")U^xU t r,\ip ). (35) 







B. Non-Hermitian coupling 

The case where the coupling operator is non-Hermitian is slightly more complicated, 
although the general argument is similar. Consider the Hamiltonian 



H tot = H + ^gj{La) + Pcij) +^2%ujja}i 

J=l 3=1 

= H + K + K ] (36) 

where K = £™ =1 g 3 La) + (1/2) £f=i hu^a,. Now suppose that we can find a complete or 
over-complete eigenbasis for L such that L\X) = A | A) and J <i 2 A|A)(A| = 1. If and \v) are 
two such states then 



(fj,\p(t)\v) = J d 2 ^' J d 2 iy'(fi'\%jjo)(i/jo\iy') J d 2 a\... J d 2 a r . 



oo co m 



e ••• e na-™>? 

ni=0 n m =0jr'=l 

(//, ai,..., aje-" 1 * **/*^', ni, ... , n m ) (j/, nt,..., n m \e iHtott/h \fi, a x , . . . , a m ) (37) 
where = e - huJ j/ k BT_ Tj s i n g Trotter product formula 

e -iH to tt/h _ e -iHdt/h e -ikdt/h e ~iK^ dt/% ^-iHdt/h ^-iKdt/h g-^ 1 " dt/h ^gg^ 

N^oo 

and inserting closure relations between the K and factors then gives 

(//, ai,..., a m |e~ iiW/a |//, ni, . . . , n m ) = lim / d 2 /ii • • • / ^Viv-i 
(p\e- i6dt / h \fi N ^) . . . (^e-^l^)^) 

m 

{i/h)[g j ^ N - 1 a\+{l/2)hu) j d\a j ]dt/h -{i/h)[gj^* N _ 1 a\ + {l/2)hu) j d\d j ]dt/h 



e 



and combining terms using Eq. (jlUJl and introducing path integral notation gives 
(fi, a^e-^^n', n 1) ...,n m )= f e 45 ^ 



The recipe for dealing with multiple coupling operators is a straightforward generalization 
of this result since each has its own oscillator bath. 
Defining ipj again via 

^j(t, a jt a*) = (ayle-^^^o *V a Wo dt'^+hu^t]^ ^ 
and differentiating then gives 

dtpj (t,aj,etj)/dt= - (if h) {gj (/^a* + /i* (atj /2+d/ da* ) ) + hujja* (a 3 -/2 + d/ da*)}i/jj (t, atj , a*) . 

(42) 

The ansatz ()19|) works here also and consistency then requires 

bj(t) = -(i/h) f dt' tig^W 
Jo 

c .(t) = -d/%) f dt' ntgje-^-^ 



o 



d 3 (t) = ~{l/h 2 ) J* dt' J* dt" nlnrgte-^-^ -irijUJjt. (43) 

Substituting this back into the subsystem density matrix element, performing the sum over 
the rij for each j, and then performing the integrals over atj gives 

W(t)\") = J rf V/ d 2 v'W )^ W) fv[pL) £v[u] e iS ^ h e- iS ^ n 

e Io dt ' So dt"n tl up,a-*(t',t")+ Jo dt' J* dt"n*v t ,,a+*(t',t") e - f* dt' // dt" f^cT (tf ,t")~ J* dt' dt" H , ^,a+ (f ,t") 



e 



So dt ' So dt»v t ,vl„ a -*(t>,t")-f Q dt' ft dt"v*v t „a+*(t',t") ( 44) 



where a -(t',t") = J2f =1 ^Jy^-e"^ '(*'-*") and a+(t',t") = Y^i^T^r/^'' 1 '^ ■ Now we 
introduce two independent complex noises and zf and unravel the coupled terms via 



e So dt ' S dt "^K" a *(*'•*") = M[eSo dt ' 0v* t , 

e So dt ' So dt"^ f u t „a+*(t',t") = M [J* dt' (n* tf z++v t ,z+*)^ ^ 

where M[ ] denotes the Gaussian average over both noises. It now follows that 

(fx\p(t)\^ = nWt)(^tW)} (46) 

where 

(^ t ) = J dWl^o) fv\ii\ e iS ^ n eSo dt ' (*vv+^*^ 

e - J* dt' f*' dt"^ t „a-{t',t")- f* dt' dt"^ tl ^„a+(t',t")_ ^ 
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This path integral representation of the wavefunction is also of dynamical semigroup form 
and hence we can introduce a propagator U t via (fi\ipt) = {fi>\Ut\ipo) which then obeys 

dU t /dt = -{i/h)HU t + LU t z; + tiu t z+ 

-L ] U t f dt"a~{t, t")U^LUt" - LU t f dt"a + (t, t")U^ti(j t „ (48) 

JO JO 

and this can again be recast in a norm preserving form. 



C. Multiple non-Hermitian couplings 

Finally, consider the case where we have p non-Hermitian coupling operators and non-zero 
temperature. Consider a Hamiltonian of the form 

H tot = H + j^{k l + k\) (49) 
1=1 

where Ki = Sj=i 9j^i^j,i + ( 1/2) Yljli where a different subset of bath operators 

couples to each system operator L\. Now assume that each Li has a complete or overcomplete 
eigenbasis Li\\i) = and / d 2 Xi\Xi) (A/| = 1. Pick eigenstates and \v) of Li, so that 

we may write 

W(t)W) = / <*V / d 2 v'(n'\^){Mv') J d 2 a\...J d 2 < p 

p oo oo mi ( 

n e e na-^Ki 

1=1^=0 n^=0i=l 

(/i, a}, . . . , al p \e-^ h \^, n{, . . . , < p ) <*/, n{, . . . , <Je^^|//, a}, . . . , < p > (50) 
where u^y = e - hu >j,ilk B T _ ^ ow we em pi y a Trotter product formula 

-ititott/h _ jj m -iHdt/h -iKidt/h -ikldt/h -iK p dt/h -iK$dt/h 

e -iHdt/h e -ikidt/fi e -ikldt/h g-ikpdt/hg-ikldt/h ^g-Q 

in which the p+ 1 factors are repeated iV times. Now we insert closure relations between each 
e -ikidt/h e -ikjdt/h p a i T This results in a sort of overcomplete path integral representation in 
which 

11 



Now, each factor can be handled as in the previous subsection. When the results are 
substituted into Eq. (|5Uj) one obtains 

[v[u l ] J V\v 2 \ ...J V[u p ] e iS ^'-^ h e -*s[ v \..,uvyn 

f[ e fo dt ' So dt"^u^-*(t',t")+f* df f* dt»^vl„*i+*(t>,t») e - dt> // dt'^ct-VW-fl dt> dt"^lW + (t',t») 

1=1 

e - So dt ' So dt>>uy t *a<-*(t>,n-f* dt' // dt»^„a'+*(t',t") (53) 

from which the stochastic wave function can be obtained by introducing two independent 
complex noises for each of the p coupling operators. The equation for the propagator can 
then be obtained along the lines followed in the previous subsections. 



III. GENERAL DYNAMICAL EQUATIONS 

Consider a subsystem-bath model with multiple operators L k interacting with different 
subsets of the radiation field 

H tot = H + 9 k j{L k a) >k + L\a j)k ) + nuJ j,kd] >k d jyk (54) 

k=ij=i fe=ij=i 

where H is the subsystem Hamiltonian, L k are system coupling operators, and #| is a 
coupling constant for an oscillator mode of frequency ojj ik of manifold k. In the norm- 
preserving formulation of NMQSD at zero temperature the evolution of the state vector ip t 
is governed by the VDE 

^ = -(i/h)H 4> t + it ~z k {L k - (L k ) t ) 
- ±fdsa\t lS ) (Li-<Lt) t )|| 

k=l J0 0Z s 

+ ± fdsa k {t,s) (MLl-(Ll) t )&)^ (55) 

k =l J0 0Z s 

where z k = z\ + J^ds a k *(t,s)(L k ) s and z\ is a complex colored noise[20] for manifold k 
with correlation function a k (t, s) = M[z k *z k ] = (l/h 2 ) Y^ =l ^e _ ^' (t_s) . The notation (L k ) t 
denotes the quantum expectation {ip t \L k \ip t ) while M[. . .] denotes an average over different 
realizations of the noises. The exact reduced density matrix p t of the subsystem is given as 
an average of diadics via p t = M[\x/j t )(ipt\]- 
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When we reformulate the theory in terms of a propagator ip t = U t ipo then we obtain a 
closed set of equations 

dfl ft 

— - = -(i/h)HU t + J2 (L k -(L k ) t )U t (z k + dsa k *(t,s)(Ll) s ) + C t U t 
dt k=1 Jo 

- ~ (Lt)t) &t [ds a k (t,s)U; l L k U s 

where 

n ft 
Ct = E^o|^(4 - {L\)t)Ut / ds a^sjU^L.U^o) (57) 
k=i Jo 

depends on the initial state ipo- 

For non-zero temperatures and non-Hermitian coupling operators the equations are 



dUt_ 
dt 



-(t/h)HU t + J2 (L k -(L k ) t )U t (zt+ I dsa k -*(t,s)(Ll) s ) 
k=i Jo 
n rt 
+ E Cl\~ (Ll)t) Ut (z k+ + / ds a k+ *(t, s)(L k ) s ) + C t U t 
k=i Jo 
n ft 

- E(4 - (Ll)t) Ut / ds a k -(t, s)U; l L k U s 
k=i Jo 

-E(4 - (L k ) t ) U t f ds a k+ (t,s)U^L{U s 
k=i Jo 



dUi 1 _ fr-idUtf,-! 



-U t -'—^U t - L (58) 



dt 1 dt 



where 



Ct = I>o|#(£J - ( L t)t)Ut fds a k -{t,s)U; l L k U s \i, ) 
k=i Jo 

+ J2(^o\U}(L k - (L k ) t )U t fds a k+ (t lS )U; l LlU s \^) (59) 
k=i Jo 

The most efficient way of solving these equations depends on the properties of the memory 

functions. We will assume that the memory functions consist of a few terms of exponential 

form, i.e., 

a k (t, s) = E A^e-^-^e-™^-^ (60) 

3=1 

where A^ k and 7^ are positive numbers. The terms do not in general correspond to 
physical bath oscillator modes. Instead the expansion can be viewed as a best fit to 
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the memory function, obtained by nonlinear least squares or other techniques. In many 
cases the number of required terms can be quite small. The case where the memory 
function cannot be represented this way is considered elsewhere p| . Defining operators 
V t j ' k = Jq(Is A jtk e-^i> k{t - s) e~ tul i' k{t - s) U- l L k U s then Eqs. ® become 

m k 



dUt 
dt 



-(i/h)HU t + £ (L k - (L k ) t ) U t (zf + £ yf ) + C t U t 

k=l j=l 
k=l j=l 



dUt 1 = ^-idEijj-i 

dt 1 dt 1 

dVi 



3,k 



dt 

ft 



dyi' k 



-{jj,k + iujj,k)V t 3,k + A jtk U t l L k U t 
df -(lj,k - i^j,k)yt k + Aj,k(Ll) t (61) 

where y{' k = J * dsAj^e' 1 ^^ 1 '^ e 1 ^^^ 1 '^ {V k ) s . Colored noises z\ = J2T=i^t' can be gener- 
ated using stochastic equations |q] 



dg' k = -(7i,fc + + ^ jtk Aj, k dWt k (62) 

which are integrated from — oo to t[5}. Here Wl' k are complex Wiener processes with proper- 
ties M[dWi' k dWi' k ] = and M[dW t jM dW l s ' m ] = 8 ttS 5 jtl 5 k>m . Accurate and efficient methods 
for solving sets of equations like (JoTj) and (jp^j) are well established 12, 3- Similar consider- 
ations apply in the case of finite temperature. 

IV. THREE-LEVEL MODEL FOR 24 Mg+ 

The six levels which comprise the relevant electronic states of 24 Mg + can be mapped to 



a three-level system 10] along the lines considered by Hulet and Wineland[8]. To do this we 
defines levels labeled 1, 2 (bright states) and 3 (dark state) from the six levels of via the 
correspondences 

Pll <— Pll + P55 

P22 *— P33 

P33 <- P22 + P44 + P66- 
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State 1 is now the ground state and 2 is an excited state lying about 4.4 eV above 1. The 1-2 

transition is resonantly driven by a laser and fluorescence from the 2 state is monitored by 

photodetectors. Occasionally interactions with the ambient radiation field cause transitions 

to state 3 which is near resonant with 1. During these periods fluorescence stops. 

1 

To model this dynamics with NMQSD we must obtain the Hulet-Wineland 8| equations 
in the Markovian limit. Near steady-state the mapping 

R-P33 *~ (2/3)7p44 
R+{pll + p22) <- (2/3)7p 55 

is valid where R- and i?+ are rates into and out of the bright manifold (1+2). These 
quantities are defined in [8( as i?_ = 8fi 2 7/9a 2 and R + = fi 2 7/18a 2 . Here Vt is the Rabi 
frequency of the laser and a is a Zeeman shift, while 7 is the spontaneous decay rate out of 
level 2. These developments together with Eqs. (1) from [8] indicate that the equations for 
the diagonal density matrix elements of the three-level system are 

p u = nimp 21 +'jp22 + R-p33~ R+(pn+ P22) 
P22 = -Olm/^i - IP22 

f>33 = -R-P33 + R+(PU + P22)- 

This is the set of equations we need to reproduce as closely as possible in the Markovian 
limit. 

A previous application of Markovian quantum state diffusion to 24 Mg + introduced a set 
of coupling operators [l^| 



u 


= A 12 |l)(2| 


(63) 


U 


= A 13 |l)(3| 


(64) 


U 


= A 31 |3)(l| 


(65) 


U 


= A(|l)(l|-|3)(3|) 


(66) 



with system Hamiltonian 

A=^(|1)<2| + |2><1|). (67) 

The parameters A12, A13, A31 and A were treated as free variables and they were chosen to fit 
experimental data. This previous study cannot therefore be considered predictive. Here we 
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will try to choose the parameters to reproduce the results of Hulet and Wineland ^] . The 
coupling operator L\ governs spontaneous emission from 2, L2 and L3 mediate transfers into 
and out of the dark state 3, and L 4 mode 
we then obtain the Lindblad-Kossakowski 



s the photo-detectors [8]. In the Markovian limit 



211 ] type equation 



^ = -(i/fi)[H,p}+Tj2[L h pLl - {l/2)L\L k p-{l/2)pL\L k } (68) 

k=l 

where we have assumed that the coupling operators share a common memory time r = 
Jq°° dt Re a(t, 0). The diagonal matrix elements then satisfy 

p u = nimp 21 + T\ 2 12 p 22 + r\ 2 13 p 33 - rA^pn 
P22 = -film/021 - T\\ 2 p 2 2 

P33 = -tA? 3 /?33 + rX^pu (69) 

which are quite similar to Eqs. (J63|) . The correspondence is not exact because Eq. (}68|) 
is positivity preserving while (fF)3*J) is not. In any case we can now compare the two sets of 
equations and deduce that 



A 12 = \h r 



A 13 = JR-/t 



A31 = V jR +/ r 

(70) 

which reduces the set of unknowns in the model. Some of the remaining constants are 
known. For example, 7 = (2tt) 43 MHz, a = (47r) 26.1 GHz for a magnetic field of 1.4 T. In 
time units of 10~ 3 7" 1 = 3.7 psec we find a = 12.1, and we set fi = 2. 

The remaining unknowns are the parameter A associated with the photodetector, and 
the temperature and distribution of frequencies for the memory function. None of these 
quantities can even be estimated. Accordingly we arbitrarily set A = .22/y / r and choose a 
memory function common to all coupling operators at K which is of the form ()6())1 with 
parameters as given in the table. This memory function is strongly non-Markovian with an 
initial fall-off at around 50 time units. The parameter r can be calculated from a(t, 0) and 
we found r = 508.6. 
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A- 


'"v ■ 




9 487407^4 


004^779QS84 


-0 0QS49^S66^ 


5 52627445 


010808938 


-0 0766453125 


10. 


0.0271137624 


0.00120546934 


9.58905445 


0.0205613891 


-0.0457549602 


8.16208273 


0.0269287619 


0.0599005113 



V. RESULTS 

We used the reformulated version of NMQSD and the Markovian QSD theory to calculate 
individual realizations of the stochastic dynamics of the driven ion. In each case we chose 
l^o) — The dynamical variable we chose to monitor was the probability of the system 
to be in the bright manifold 

p(0HW t )l 2 + IW t >l 2 (7i) 

which is analogous to the experimentally observed fluorescence intensity. A few example 
trajectories are shown in Fig. 1. Note that the dark periods tend to be more frequent and 
shorter for the Markovian dynamics. The Markovian trajectories are also more noisy as one 
would expect. 

It is customary in such jump experiments to construct a histogram of the frequency with 
which fluorescence intensities are observed. Mathematically the histogram x(P) vs P is 
given as a limit 

X (P) = A l|m o ((l/T) J o dt J p dy 5(y - P(t))) (72) 

which we approximate by choosing a finite but small AP. The angle brackets denote an 
ensemble average over individual trajectories (4000 in the non-Markovian case and 10000 in 
the Markovian case). For 24 Mg + this histogram was observed to have two peaks, one near 
zero signal strength and one near maximum signal strength. The ratio of the area under 
the maximum signal peak to that under the minimum signal peak is 16 according to both 
theory and experiment. The numerical results we obtained are shown in Fig. 2. In the 
non-Markovian case two clear peaks are observed, one near P = and one near P — 1, 
corresponding to occupation of the dark and bright manifolds respectively. To extract a 
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(a) 



0.015625 
0.00390625 
0.000976562 
0.000244141 




(b) 

FIG. 2: Non-Markovian (a) and Markovian (b) signal frequency x(-P) vs P 
ratio of areas under the peaks we fit the data to 

X (P) = Xl(P) + X2(P) + Xbackground(P) 

where 

Xi(P) = h 1 VP/{(P/w 1 ) 2 + l} 
models the low signal peak and 

X2(P) = h 2 P/{[(l - P)/w 2 f + [(I - P)/w 3 ] 2 + (I - P)/w A + 1} 
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+h 3 P/{[(l - P)/w 5 } 2 + (1 - P)/w 6 + 1} 



(75) 



models the high signal peak. The background in between was fitted to 

Xback 9 round(P) = h 3 {eii[(P - P 1 )/w 7 ] - erf[(P - P 2 )/w 8 ]}. (76) 

Using this fitting function we were able to extract a ratio of 16.1 for the area under the 
strong signal peak to the area under the weak signal peak. The best fit is shown in Fig. 3. 

In the Markovian case the existence of the low signal peak is less clear. We chose fitting 
functions of the form 

X i(P) = h 1 P/{(P/w 1 ) 2 + P/w 2 + 1} (77) 

for the low signal peak and 

X2(P) = h 2 P/{[{\ - P)/w 3 ] 2 + (1 - P)/w A + 1} 

+h 3 P/{[(l - P)/w 5 ] 2 + (1 - P)/w 6 + 1} (78) 

for the high signal peak, which are similar to those of the non-Markovian case. The back- 
ground model function was the same as in the non-Markovian case. The ratio of areas 
extracted was 16.8. This minor discrepancy is likely due to lack of Monte-Carlo conver- 
gence. 

Thus, we are able to verify the ratio of areas in both the Markovian and non-Markovian 
cases. The mathematical form of the histogram seems to vary little from the Markovian to 
non-Markovian cases even though the parameters and appearance of the two distributions 
are quite different. This is encouraging since it may be possible in future to obtain time 
resolved experimental histograms which could be compared to theory to see whether there 
is agreement of mathematical forms. Currently, lack of time resolution and poor detection 
efficiency in the experiments prevent detailed comparisons. 

VI. SUMMARY 

We have shown in considerable detail that NMQSD can be reformulated in terms of a 
solvable stochastic integrodifferential equation. The reformulated theory has been tested 
against exact results for a number of problems, and employed to investigate a number of 
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0.025 
P 



(a) 




(b) 

FIG. 3: Non-Markovian signal frequency x(-P) vs P at small and large signal values (solid curve) 
and best fit (dashed curve) 

problems for which exact solutions are not known &\. We expect that the method will prove 
useful for many few-level-system problems in quantum optics. 

In this manuscript we applied the theory to the problem of intermittent fluorescence 
in 24 Mg + . Previous applications of Markovian quantum state diffusion to quantum jumps 
in 24 Mg + treated all parameters as free variables (lo|, and hence cannot be considered pre- 
dictive. Our more careful study shows that quantum jumps do indeed occur on picosecond 
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timescales as had been speculated earlier |lCfl. We also computed the probability distribution 
function for signal intensity and obtained the characteristic two peaks corresponding to flu- 
orescence on and off. Our results show that both Markovian and non-Markovian versions of 
NMQSD reproduce the experimental and theoretical result that the ratio of the area under 
the bright peak is 16 times that of the area under the dark peak. We found that while the 
Markovian and non-Markovian distribution functions looked qualitatively different, they in 
fact share a very similar mathematical form. This raises a number of interesting possibilities. 
First, if time resolved experiments are possible it may be possible to verify the shape the 
histogram experimentally. This would be a much stronger result than the simple ratio of 
areas. Secondly, it raises the possibility of using the jump experiment to measure properties 
of the radiation field, since the qualitative shape of the histogram is sensitive to the memory 
function which contains quite a lot of information about the distribution of frequencies and 
temperature. 

The authors acknowledge the support of the Natural Sciences and Engineering Research 
Council of Canada. 
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